# PROLOG ##############################################################################
# PROJECT: Priority Effects Experiment
# PURPOSE: Determine how species arrival via seed additions affects the community assembly of restored prairie plant communities
#
# AUTHOR: Katherine Wynne (Wynnekat@msu.edu)
# COLLAB: L. Sullivan
# CREATED: 06 December 2022
#
# FILES: 1) inner_cover_summer_2022.xlsx - vegetative cover taken in the Summer (end of June)
# 2) inner_cover_fall _2022.xlsx - vegetative cover taken at peak biomass (end ofAugust)
# 3) Checklist of Missouri Flora.xlsx - contains the information about all plants found in MO
# 4) Cleaned_PE_Detailed_Species_List_2022.xlsx - dataset containing the information about only plants found in the study
#
#LAST UPDATED: 06 September 2023
# PROLOG ##############################################################################
NON = None (no species seeded)
SIM = Simulataneous (36 species seeded on March 22nd, 2021)
LE = Lumped early (18 early-dispersing species seeded on March 23rd, 2021 then 18 late-dispersing species seeded on September 5th, 2021)
LL = Lumped late (18 late-dispersing species seeded on March 23rd, 2021 then 18 early-dispersing species seeded on September 5th, 2021)
NAT = Natural; species added according to peak dispersal date
- May 24th 2021 (VIOPED, PACPLA)
- June 4th 2021 (SISCAM, SPHOBT, CORLAN)
- June 22nd 2021 (LOBSPI, TRAOHI, CARBUS, HEURIC)
- July 11th 2021 (CORPAL, DODMEA)
- July 25th 2021 (KOEMAC, AMOCAN)
- August 2nd 2021 (LINSUL, MELVIR)
- August 21st 2021 (DALCAN, DALPUR, ACHMIL)
- September 17th 2021 (CROSAG)
- October 4th 2021 (CHAFAS, MONFIS, RATPIN, SPOHET, RUDHIR)
- October 18th 2021 (BIDARI, HELMOL, LESCAP)
- November 1st 2021 (PENDIG, ERYYUC, HYPPUN, SORNUT, SCHSCO, LIAPYC)
- November 19th 2021 (CORTRI, PYCTEN, SOLRIG)
### Remove unnecessary columns
#### Get rid of the notes, unknown, and senensced columns
### Fall
inner_cover_fall_2021 <- inner_cover_fall_2021[,-c(6,7,8)]
#### change log
# - senesced PLASPP changed to PLAVIR since PLALAN and PLAMAJ do not senesce until Oct. (Sep 6, 2023)
# - JUNSPP changed to JUNINT since JUNINT has been observed in 2022 and 2023 (Sep 6, 2023)
### Get rid of EMP
### Fall
inner_cover_fall_2021 <- subset(inner_cover_fall_2021, Treatment != "EMP")
# Make a unique identifier for year and season
### 2021
#### Season
inner_cover_fall_2021$Season <- rep("Fall", nrow=(inner_cover_fall_2021))
#### Year
inner_cover_fall_2021$Year <- rep("2021", nrow=(inner_cover_fall_2021))
### 2022
# Already has the identifiers
inner_cover_fall_2022$Year <- as.factor(inner_cover_fall_2022$Year)
inner_cover_summer_2022$Year <- as.factor(inner_cover_summer_2022$Year)
### 2023
#### Season
inner_cover_summer_2023$Season <- rep("Summer", nrow=(inner_cover_summer_2023))
inner_cover_fall_2023$Season <- rep("Fall", nrow=(inner_cover_fall_2023))
#### Year
inner_cover_fall_2023$Year <- rep("2023", nrow=(inner_cover_fall_2023))
inner_cover_fall_2023$Year <- as.factor(inner_cover_fall_2023$Year)
inner_cover_summer_2023$Year <- rep("2023", nrow=(inner_cover_summer_2023))
inner_cover_summer_2023$Year <- as.factor(inner_cover_summer_2023$Year)
# Join years together
full_2021_data <- inner_cover_fall_2021
full_2022_data <- full_join(inner_cover_summer_2022,inner_cover_fall_2022)
## Joining with `by = join_by(Season, Year, Block, Treatment, SPP6, Percent_Cover,
## Scientific_Name)`
full_2023_data <- full_join(inner_cover_summer_2023,inner_cover_fall_2023)
## Joining with `by = join_by(Block, Treatment, SPP6, Percent_Cover,
## Scientific_Name, Season, Year)`
# Join all years together
### First 2021 and 2022
full_21_22_data <- full_join(full_2021_data, full_2022_data)
## Joining with `by = join_by(Block, Treatment, SPP6, Percent_Cover,
## Scientific_Name, Season, Year)`
### Then 2021 and 2022 with 2023
full_year_data <- full_join(full_21_22_data, full_2023_data)
## Joining with `by = join_by(Block, Treatment, SPP6, Percent_Cover,
## Scientific_Name, Season, Year)`
# lump certain taxa together
### Lump the Carex for now (*****Ask Lauren about what to do about when the CARSPP > CARFRA or CARBRE?***)
## MELOFF + MELALB -> MELSPP
## LEPVIR -> LEPSPP
## CARFRA -> CARSPP
## CARBRE -> CARSPP
for(i in 1:nrow(full_year_data)) {
if(full_year_data[i,3] == "CARFRA"){full_year_data [i,3] <- "CARSPP"}
if(full_year_data[i,3] == "CARBRE"){full_year_data [i,3] <- "CARSPP"}
if(full_year_data[i,3] == "MELOFF"){full_year_data [i,3] <- "MELSPP"}
if(full_year_data[i,3] == "MELALB"){full_year_data [i,3] <- "MELSPP"}
if(full_year_data[i,3] == "LEPVIR"){full_year_data [i,3] <- "LEPSPP"}
# Make all cover data that was less than 1 = to 1 (to make the data discrete)
if(full_year_data[i,4] < 1) {full_year_data [i,4] <- 1}
}
inner_cover_max <- full_year_data %>%
group_by(Block, Treatment, Year, SPP6) %>%
summarise(max_cover=max(Percent_Cover))
## `summarise()` has grouped output by 'Block', 'Treatment', 'Year'. You can
## override using the `.groups` argument.
##### Note I manually checked species that were found in both the datasets (ACHMIL, TRIREP, CORLAN) and confirmed that the code above took the highest cover value for a species seen twice
### Remove Bare
inner_cover_max_only <- subset(inner_cover_max, SPP6 != "BARE")
### Remove Litter
inner_cover_max_only <- subset(inner_cover_max_only, SPP6 != "LITTER")
head(inner_cover_max_only)
## # A tibble: 6 × 5
## # Groups: Block, Treatment, Year [1]
## Block Treatment Year SPP6 max_cover
## <dbl> <chr> <chr> <chr> <dbl>
## 1 1 LE 2021 ACAVIR 1
## 2 1 LE 2021 ACHMIL 1
## 3 1 LE 2021 AMATUB 1
## 4 1 LE 2021 AMBART 20
## 5 1 LE 2021 BARVUL 1
## 6 1 LE 2021 BROJAP 1
### Bare dataset
inner_bare_cover <- subset(inner_cover_max, SPP6 == "BARE")
### Litter dataset
inner_litter_cover <- subset(inner_cover_max, SPP6 == "LITTER")
### Manipulate 2021 dataset to match 2022 and 2023
#### Bare Ground
bare_2021 <- subset(bare_cover_2021, Cover_Type != "Litter")
bare_2021 <- subset(bare_2021, Treatment != "EMP")
names(bare_2021) <- c("Block", "Treatment", "SPP6", "max_cover")
bare_2021$Year <- rep("2021", nrow(bare_2021))
for(i in 1:nrow(bare_2021)) {
if(bare_2021[i,3] == "Bare"){bare_2021[i,3] <- "BARE"}
}
#### Litter
litter_2021 <- subset(bare_cover_2021 , Cover_Type != "Bare")
litter_2021 <- subset(litter_2021, Treatment != "EMP")
names(litter_2021) <- c("Block", "Treatment", "SPP6", "max_cover")
litter_2021$Year <- rep("2021", nrow(litter_2021))
for(i in 1:nrow(litter_2021)) {
if(bare_2021[i,3] == "Litter"){bare_2021[i,3] <- "LITTER"}
}
### Full_join datasets
inner_bare_cover <- full_join(inner_bare_cover, bare_2021)
## Joining with `by = join_by(Block, Treatment, Year, SPP6, max_cover)`
inner_litter_cover <- full_join(inner_litter_cover, litter_2021)
## Joining with `by = join_by(Block, Treatment, Year, SPP6, max_cover)`
### Add species information to the dataset
#### N/A is N = native, A = non-native, G = Genus
inner_cover_max_only <- full_join(inner_cover_max_only, species_list_PE)
## Joining with `by = join_by(SPP6)`
#Remove unnecessary columns
inner_cover_max_reduced <- inner_cover_max_only[, -c(6,7,10)]
Below is a summary table showing the mean and SD for mean C value (C), species richness (SR), and Shannon diversity index (SDI) for each seeding treatment.
Distribution of C values for all treatments were right skewed, indicating the majority of native plants in all treatments were of low conservatism value. SIM had the most “even” distribution of C values and almost every value was represented in this treatment.
## Warning in geom_density(bins = 30, alpha = 0.4): Ignoring unknown parameters:
## `bins`
# A couple of outliers in the NAT treatment but likely alright since the variablity in outcomes seems to be a characteristic of NAT
# only one extreme value in LE 2021 (likely because of DALCAN and VIOPED)
Summary_table %>%
group_by(Treatment, Year) %>%
identify_outliers(Mean_C)
## # A tibble: 4 × 8
## Treatment Year Block Species_richness Shannon Mean_C is.outlier is.extreme
## <chr> <chr> <dbl> <int> <dbl> <dbl> <lgl> <lgl>
## 1 NAT 2023 4 18 2.44 2.69 TRUE FALSE
## 2 NON 2023 6 19 1.99 1.18 TRUE FALSE
## 3 SIM 2022 5 20 1.83 3.13 TRUE FALSE
## 4 SIM 2023 2 27 2.41 3.42 TRUE FALSE
# Mean C fits a normal distribution
Summary_table %>%
group_by(Treatment, Year) %>%
shapiro_test(Mean_C)
## # A tibble: 10 × 5
## Treatment Year variable statistic p
## <chr> <chr> <chr> <dbl> <dbl>
## 1 LE 2022 Mean_C 0.926 0.548
## 2 LE 2023 Mean_C 0.972 0.907
## 3 LL 2022 Mean_C 0.890 0.317
## 4 LL 2023 Mean_C 0.861 0.193
## 5 NAT 2022 Mean_C 0.901 0.379
## 6 NAT 2023 Mean_C 0.935 0.616
## 7 NON 2022 Mean_C 0.928 0.566
## 8 NON 2023 Mean_C 0.882 0.280
## 9 SIM 2022 Mean_C 0.816 0.0809
## 10 SIM 2023 Mean_C 0.879 0.263
ggplot(data = Summary_table , aes(x= Mean_C)) +geom_histogram(bins =15)
# Vast majority of points fall along the reference line
ggqqplot(Summary_table , "Mean_C", ggtheme = theme_bw()) +
facet_grid(Year~ Treatment, labeller = "label_both")
Mean C values were significantly higher in all treatments compared to the NON treatment.
#Dropped the mixed model due to singular fit -> the variance of the random effect block is estimated close to zero and does not further inform the data.
Summary_table$Block <- as.factor(Summary_table$Block)
# Normal interaction model
#Mean_C.mod.lmer.interaction <- lmer(Mean_C~Treatment*Year+(1|Block), data =Summary_table)
#AIC(Mean_C.mod.lmer.interaction)
# Normal - additive model (has the lowest AIC)
Mean_C.mod.lm.additive <- lm(Mean_C~Treatment+Year, data =Summary_table)
## Summary results
summary(Mean_C.mod.lm.additive)
##
## Call:
## lm(formula = Mean_C ~ Treatment + Year, data = Summary_table)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.99511 -0.34068 -0.07964 0.27121 1.09720
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.6743 0.1583 10.575 9.11e-15 ***
## TreatmentLL 0.4931 0.2044 2.412 0.019270 *
## TreatmentNAT -0.3604 0.2044 -1.763 0.083548 .
## TreatmentNON -1.3663 0.2044 -6.684 1.34e-08 ***
## TreatmentSIM 0.7765 0.2044 3.799 0.000371 ***
## Year2023 0.2812 0.1293 2.175 0.034032 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5007 on 54 degrees of freedom
## Multiple R-squared: 0.7199, Adjusted R-squared: 0.694
## F-statistic: 27.76 on 5 and 54 DF, p-value: 8.389e-14
## AIC test
AIC(Mean_C.mod.lm.additive)
## [1] 94.93578
## ANOVA type 3
anova(Mean_C.mod.lm.additive)
## Analysis of Variance Table
##
## Response: Mean_C
## Df Sum Sq Mean Sq F value Pr(>F)
## Treatment 4 33.612 8.4029 33.5209 4.716e-14 ***
## Year 1 1.186 1.1858 4.7304 0.03403 *
## Residuals 54 13.537 0.2507
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## I removed the random effect block because it had no effect on the model fit and explained no additional variance. (chose the most parsimonius model)
# r.squaredGLMM(Mean_C.mod.lmer.additive)
## Inclusion of the random effect marginally increased model fit
est.lm.treatment <- emmeans::emmeans(Mean_C.mod.lm.additive, "Treatment")
#Posthoc test
## Treatment
pairs(est.lm.treatment, adjust = "tukey")
## contrast estimate SE df t.ratio p.value
## LE - LL -0.493 0.204 54 -2.412 0.1276
## LE - NAT 0.360 0.204 54 1.763 0.4054
## LE - NON 1.366 0.204 54 6.684 <.0001
## LE - SIM -0.777 0.204 54 -3.799 0.0033
## LL - NAT 0.853 0.204 54 4.176 0.0010
## LL - NON 1.859 0.204 54 9.097 <.0001
## LL - SIM -0.283 0.204 54 -1.387 0.6388
## NAT - NON 1.006 0.204 54 4.921 0.0001
## NAT - SIM -1.137 0.204 54 -5.562 <.0001
## NON - SIM -2.143 0.204 54 -10.483 <.0001
##
## Results are averaged over the levels of: Year
## P value adjustment: tukey method for comparing a family of 5 estimates
# Residual plot looks decent
plot(Mean_C.mod.lm.additive)
Summary_table_mean$Treatment <- factor(Summary_table_mean$Treatment, levels=c("NON", "NAT", "LE", "LL", "SIM"))
lower.c <- Summary_table_mean$C_Mean - Summary_table_mean$C_SD
upper.c <- Summary_table_mean$C_Mean + Summary_table_mean$C_SD
pd <- position_dodge(width = 0.8)
C.plot <- ggplot(data = Summary_table_mean, aes(x = Year, y = C_Mean, color = Treatment))+
geom_point(data = Summary_table_mean, size = 3, aes(shape = Treatment, fill = Treatment), position = pd)+
geom_errorbar(aes(ymin = lower.c , ymax = upper.c), width = 0.5, position = pd) +
labs(y = "Mean C", x = "Year" )+
theme_classic() +
theme(text=element_text(size=16), legend.key.size=unit(1, "cm"), legend.position="none")+
scale_color_manual(name = "Treatment", values=c("#E69F00", "#F0E442", "#009E73", "#56B4E9", "#0072B2"), breaks=c("NON", "NAT", "LE", "LL", "SIM"), labels = c("NON", "NAT", "LE", "LL", "SIM"))+
scale_fill_manual(name = "Treatment", values=c("#E69F00", "#F0E442", "#009E73", "#56B4E9", "#0072B2"), breaks=c("NON", "NAT", "LE", "LL", "SIM"), labels = c("NON", "NAT", "LE", "LL", "SIM"))+
scale_shape_manual(name = "Treatment", values=c(21, 22,23, 24, 25), breaks=c("NON", "NAT", "LE", "LL", "SIM"), labels = c("NON", "NAT", "LE", "LL", "SIM"))
C.plot
SIM treatments had the greatest species richness compared to all other treatments, while NON and NAT had the lowest.
# A couple of outliers in the NAT treatment but likely alright since the variablity in outcomes seems to be a characteristic of NAT
Summary_table %>%
group_by(Treatment, Year) %>%
identify_outliers(Species_richness)
## # A tibble: 5 × 8
## Treatment Year Block Species_richness Shannon Mean_C is.outlier is.extreme
## <chr> <chr> <fct> <int> <dbl> <dbl> <lgl> <lgl>
## 1 LE 2023 3 15 2.10 2.2 TRUE FALSE
## 2 LL 2022 2 25 2.31 2.06 TRUE FALSE
## 3 NAT 2022 4 22 2.18 2.31 TRUE TRUE
## 4 NAT 2022 6 15 2.10 1.09 TRUE FALSE
## 5 NON 2023 6 19 1.99 1.18 TRUE FALSE
# Species richness fits a normal distribution
Summary_table %>%
group_by(Treatment, Year) %>%
shapiro_test(Species_richness)
## # A tibble: 10 × 5
## Treatment Year variable statistic p
## <chr> <chr> <chr> <dbl> <dbl>
## 1 LE 2022 Species_richness 0.996 0.998
## 2 LE 2023 Species_richness 0.847 0.149
## 3 LL 2022 Species_richness 0.814 0.0778
## 4 LL 2023 Species_richness 0.825 0.0969
## 5 NAT 2022 Species_richness 0.836 0.122
## 6 NAT 2023 Species_richness 0.884 0.286
## 7 NON 2022 Species_richness 0.960 0.817
## 8 NON 2023 Species_richness 0.955 0.784
## 9 SIM 2022 Species_richness 0.893 0.332
## 10 SIM 2023 Species_richness 0.873 0.238
ggplot(data = Summary_table, aes(x= Species_richness)) +geom_histogram(bins =15)
# Vast majority of points fall along the reference line
ggqqplot(Summary_table, "Species_richness", ggtheme = theme_bw()) +
facet_grid(Year~ Treatment, labeller = "label_both")
# Model Fitting
## Despite being count data, the distribution of Species richness is remarkabley balanced and not skewed. Model comparison using AIC reveals that the Normal model works better than the Poisson model. Model results also differ likely because of under dispersion. The more complicated Quasipoisson model that accounts for under dispersion produces similar results to the Normal mixed model. Therefore, I'm going to go with the less complicated Normal model.
SR.mod.normal <- lmer(Species_richness~Treatment+Year+(1|Block), data =Summary_table)
summary(SR.mod.normal)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Species_richness ~ Treatment + Year + (1 | Block)
## Data: Summary_table
##
## REML criterion at convergence: 283.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3018 -0.7002 0.0875 0.6078 1.6946
##
## Random effects:
## Groups Name Variance Std.Dev.
## Block (Intercept) 2.378 1.542
## Residual 7.355 2.712
## Number of obs: 60, groups: Block, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 19.9333 1.0639 20.8212 18.736 1.64e-14 ***
## TreatmentLL 0.3333 1.1072 49.0000 0.301 0.764634
## TreatmentNAT -1.7500 1.1072 49.0000 -1.581 0.120398
## TreatmentNON -4.3333 1.1072 49.0000 -3.914 0.000280 ***
## TreatmentSIM 4.2500 1.1072 49.0000 3.839 0.000355 ***
## Year2023 -1.7000 0.7002 49.0000 -2.428 0.018911 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TrtmLL TrtNAT TrtNON TrtSIM
## TreatmentLL -0.520
## TreatmntNAT -0.520 0.500
## TreatmntNON -0.520 0.500 0.500
## TreatmntSIM -0.520 0.500 0.500 0.500
## Year2023 -0.329 0.000 0.000 0.000 0.000
## Treatment and Year significant
## Interaction between Treatment and Year not significant
anova(SR.mod.normal, type = 3)
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Treatment 474.77 118.69 4 49 16.1381 1.704e-08 ***
## Year 43.35 43.35 1 49 5.8941 0.01891 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Marginal R^2 (variance explained by ONLY fixed effects) = 0.4547589
## Conditional R^2 (variance explained by the entire model) = 0.5554578
## Inclusion of the random effect increased model fit
r.squaredGLMM(SR.mod.normal)
## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help page.
## R2m R2c
## [1,] 0.4743091 0.6027586
#Posthoc test
est.lmer.SR.treatment <- emmeans::emmeans(SR.mod.normal, "Treatment")
## Treatment
pairs(est.lmer.SR.treatment , adjust = "tukey")
## contrast estimate SE df t.ratio p.value
## LE - LL -0.333 1.11 49 -0.301 0.9981
## LE - NAT 1.750 1.11 49 1.581 0.5166
## LE - NON 4.333 1.11 49 3.914 0.0025
## LE - SIM -4.250 1.11 49 -3.839 0.0031
## LL - NAT 2.083 1.11 49 1.882 0.3406
## LL - NON 4.667 1.11 49 4.215 0.0010
## LL - SIM -3.917 1.11 49 -3.538 0.0076
## NAT - NON 2.583 1.11 49 2.333 0.1519
## NAT - SIM -6.000 1.11 49 -5.419 <.0001
## NON - SIM -8.583 1.11 49 -7.753 <.0001
##
## Results are averaged over the levels of: Year
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 5 estimates
# Poisson Model
# SR.mod.pois <- glmer(Species_richness~Treatment+Year+(1|Block),family = "poisson", data =Summary_table)
# summary(SR.mod.pois)
# Anova(SR.mod.pois, type = 3)
# Check dispersion
## Underdispersed! 0.659
## Interaction between treatment and year are not significant
# Dispersion_glmer(SR.mod.pois)
# Not accounting for the underdispersion mattered a lot! .
## Used a quassipoisson model instead and year is significant
#m6 <-glmmPQL(Species_richness~Treatment+Year,
# random = ~ 1 | Block,
# family = quasipoisson(link='log'),
# data = Summary_table)
# summary(m6)
# Anova(m6, type = 3)
# Residual plot looks decent except for one large outlier
plot(SR.mod.normal)
Summary_table_mean$Treatment <- factor(Summary_table_mean$Treatment, levels=c("NON", "NAT", "LE", "LL", "SIM"))
lower.SR.c <- Summary_table_mean$SR_Mean - Summary_table_mean$SR_SD
upper.SR.c <- Summary_table_mean$SR_Mean + Summary_table_mean$SR_SD
pd <- position_dodge(width = 0.8)
SR.plot1 <- ggplot(data = Summary_table_mean, aes(x = Year, y = SR_Mean, color = Treatment))+
geom_point(data = Summary_table_mean, size = 3, aes(shape = Treatment, fill = Treatment), position = pd)+
geom_errorbar(aes(ymin = lower.SR.c , ymax = upper.SR.c), width = 0.5, position = pd) +
labs(y = "Total Species Richness", x = "Year" )+
theme_classic() +
ylim(10,30)+
theme(text=element_text(size=16), legend.key.size=unit(1, "cm"))+
scale_color_manual(name = "Treatment", values=c("#E69F00", "#F0E442", "#009E73", "#56B4E9", "#0072B2"), breaks=c("NON", "NAT", "LE", "LL", "SIM"), labels = c("NON", "NAT", "LE", "LL", "SIM"))+
scale_fill_manual(name = "Treatment", values=c("#E69F00", "#F0E442", "#009E73", "#56B4E9", "#0072B2"), breaks=c("NON", "NAT", "LE", "LL", "SIM"), labels = c("NON", "NAT", "LE", "LL", "SIM"))+
scale_shape_manual(name = "Treatment", values=c(21, 22,23, 24, 25), breaks=c("NON", "NAT", "LE", "LL", "SIM"), labels = c("NON", "NAT", "LE", "LL", "SIM"))+
theme(legend.position="bottom")
SR.plot <- SR.plot1 +
theme(legend.position="none")
leg <- get_legend(SR.plot1)
SIM and LL treatments had the greatest Shannon diversity index value while NON had the lowest. NAT and LE appear comparable.
# Two extreme outliers (one in NAT 2021 and another in SIM 2021)
Summary_table %>%
group_by(Treatment, Year) %>%
identify_outliers(Shannon)
## # A tibble: 1 × 8
## Treatment Year Block Species_richness Shannon Mean_C is.outlier is.extreme
## <chr> <chr> <fct> <int> <dbl> <dbl> <lgl> <lgl>
## 1 NON 2022 3 11 1.55 0.429 TRUE FALSE
# Except for one point, Shannon diversity fits a normal distribution
Summary_table %>%
group_by(Treatment, Year) %>%
shapiro_test(Shannon)
## # A tibble: 10 × 5
## Treatment Year variable statistic p
## <chr> <chr> <chr> <dbl> <dbl>
## 1 LE 2022 Shannon 0.971 0.901
## 2 LE 2023 Shannon 0.967 0.873
## 3 LL 2022 Shannon 0.977 0.936
## 4 LL 2023 Shannon 0.854 0.169
## 5 NAT 2022 Shannon 0.950 0.743
## 6 NAT 2023 Shannon 0.931 0.586
## 7 NON 2022 Shannon 0.848 0.151
## 8 NON 2023 Shannon 0.950 0.739
## 9 SIM 2022 Shannon 0.872 0.235
## 10 SIM 2023 Shannon 0.940 0.663
ggplot(data = Summary_table, aes(x= Shannon)) +geom_histogram(bins =15)
# Vast majority of points fall along the reference line
ggqqplot(Summary_table, "Shannon", ggtheme = theme_bw()) +
facet_grid(Year~ Treatment, labeller = "label_both")
Shannon.lmer.mod <- lmer(Shannon~Treatment+Year+(1|Block), data =Summary_table)
summary(Shannon.lmer.mod)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Shannon ~ Treatment + Year + (1 | Block)
## Data: Summary_table
##
## REML criterion at convergence: 2.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.60837 -0.47607 0.08046 0.51127 2.04855
##
## Random effects:
## Groups Name Variance Std.Dev.
## Block (Intercept) 0.007995 0.08942
## Residual 0.042202 0.20543
## Number of obs: 60, groups: Block, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.16123 0.07452 28.51199 29.003 < 2e-16 ***
## TreatmentLL 0.05159 0.08387 49.00000 0.615 0.54133
## TreatmentNAT -0.12219 0.08387 49.00000 -1.457 0.15152
## TreatmentNON -0.34644 0.08387 49.00000 -4.131 0.00014 ***
## TreatmentSIM 0.12917 0.08387 49.00000 1.540 0.12994
## Year2023 -0.01536 0.05304 49.00000 -0.290 0.77333
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TrtmLL TrtNAT TrtNON TrtSIM
## TreatmentLL -0.563
## TreatmntNAT -0.563 0.500
## TreatmntNON -0.563 0.500 0.500
## TreatmntSIM -0.563 0.500 0.500 0.500
## Year2023 -0.356 0.000 0.000 0.000 0.000
anova(Shannon.lmer.mod, type = 3)
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Treatment 1.65266 0.41316 4 49 9.7902 6.691e-06 ***
## Year 0.00354 0.00354 1 49 0.0839 0.7733
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# No significant interaction between treatment and year
## Marginal R^2 (variance explained by ONLY fixed effects) = 0.330704
## Conditional R^2 (variance explained by the entire model) = 0.3844412
## Inclusion of the random effect increased model fit
r.squaredGLMM(Shannon.lmer.mod)
## R2m R2c
## [1,] 0.358652 0.4608049
# Post-Hoc comparison
pairs(emmeans::emmeans(Shannon.lmer.mod, "Treatment"))
## contrast estimate SE df t.ratio p.value
## LE - LL -0.0516 0.0839 49 -0.615 0.9720
## LE - NAT 0.1222 0.0839 49 1.457 0.5947
## LE - NON 0.3464 0.0839 49 4.131 0.0013
## LE - SIM -0.1292 0.0839 49 -1.540 0.5420
## LL - NAT 0.1738 0.0839 49 2.072 0.2486
## LL - NON 0.3980 0.0839 49 4.746 0.0002
## LL - SIM -0.0776 0.0839 49 -0.925 0.8858
## NAT - NON 0.2242 0.0839 49 2.674 0.0727
## NAT - SIM -0.2514 0.0839 49 -2.997 0.0331
## NON - SIM -0.4756 0.0839 49 -5.671 <.0001
##
## Results are averaged over the levels of: Year
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 5 estimates
# A little heteroscedastic (bowtie shaped but probably ok)
plot(Shannon.lmer.mod)
Summary_table_mean$Treatment <- factor(Summary_table_mean$Treatment, levels=c("NON", "NAT", "LE", "LL", "SIM"))
lower.SDI.c <- Summary_table_mean$SDI_Mean - Summary_table_mean$SDI_SD
upper.SDI.c <- Summary_table_mean$SDI_Mean + Summary_table_mean$SDI_SD
pd <- position_dodge(width = 0.8)
SDI.plot <- ggplot(data = Summary_table_mean, aes(x = Year, y = SDI_Mean, color = Treatment))+
geom_point(data =Summary_table_mean, size = 3, aes(shape = Treatment, fill = Treatment), position = pd)+
geom_errorbar(aes(ymin = lower.SDI.c , ymax = upper.SDI.c), width = 0.5, position = pd) +
labs(y = "SDI", x = "Year" )+
theme_classic() +
theme(text=element_text(size=16), legend.key.size=unit(1, "cm"), legend.position="none")+
scale_color_manual(name = "Treatment", values=c("#E69F00", "#F0E442", "#009E73", "#56B4E9", "#0072B2"), breaks=c("NON", "NAT", "LE", "LL", "SIM"), labels = c("NON", "NAT", "LE", "LL", "SIM"))+
scale_fill_manual(name = "Treatment", values=c("#E69F00", "#F0E442", "#009E73", "#56B4E9", "#0072B2"), breaks=c("NON", "NAT", "LE", "LL", "SIM"), labels = c("NON", "NAT", "LE", "LL", "SIM"))+
scale_shape_manual(name = "Treatment", values=c(21, 22,23, 24, 25), breaks=c("NON", "NAT", "LE", "LL", "SIM"), labels = c("NON", "NAT", "LE", "LL", "SIM"))
SDI.plot
note some of the species used in the seed mix were present in the seed bank at my study site (e.g., RUDHIR and PENDIG). Additionally, some spillover occurred with seeded species seeding into other treatments like NON (e.g., BIDARI)
##
## Call:
## glm(formula = cbind(seeded_cover, tot_cover) ~ Treatment * Year,
## family = binomial, data = Total_Cover_Tot)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -6.4760 -2.0224 -0.4559 1.5998 4.5960
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.11998 0.08881 -23.870 < 2e-16 ***
## TreatmentLL 0.94324 0.10412 9.059 < 2e-16 ***
## TreatmentNAT -0.38112 0.13799 -2.762 0.00575 **
## TreatmentSIM 1.06270 0.10397 10.221 < 2e-16 ***
## Year2023 0.62783 0.11282 5.565 2.63e-08 ***
## TreatmentLL:Year2023 -0.41302 0.13685 -3.018 0.00254 **
## TreatmentNAT:Year2023 -0.36533 0.18174 -2.010 0.04441 *
## TreatmentSIM:Year2023 -0.37777 0.13662 -2.765 0.00569 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 837.79 on 47 degrees of freedom
## Residual deviance: 293.23 on 40 degrees of freedom
## AIC: 557.58
##
## Number of Fisher Scoring iterations: 5
## Analysis of Deviance Table (Type III tests)
##
## Response: cbind(seeded_cover, tot_cover)
## LR Chisq Df Pr(>Chisq)
## Treatment 272.445 3 < 2.2e-16 ***
## Year 32.035 1 1.514e-08 ***
## Treatment:Year 10.305 3 0.01614 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## contrast odds.ratio SE df null z.ratio p.value
## LE Year2022 / LL Year2022 0.389 0.0405 Inf 1 -9.059 <.0001
## LE Year2022 / NAT Year2022 1.464 0.2020 Inf 1 2.762 0.1048
## LE Year2022 / SIM Year2022 0.346 0.0359 Inf 1 -10.221 <.0001
## LE Year2022 / LE Year2023 0.534 0.0602 Inf 1 -5.565 <.0001
## LE Year2022 / LL Year2023 0.314 0.0328 Inf 1 -11.076 <.0001
## LE Year2022 / NAT Year2023 1.126 0.1469 Inf 1 0.909 0.9853
## LE Year2022 / SIM Year2023 0.269 0.0281 Inf 1 -12.574 <.0001
## LL Year2022 / NAT Year2022 3.760 0.4466 Inf 1 11.150 <.0001
## LL Year2022 / SIM Year2022 0.887 0.0680 Inf 1 -1.558 0.7751
## LL Year2022 / LE Year2023 1.371 0.1210 Inf 1 3.572 0.0085
## LL Year2022 / LL Year2023 0.807 0.0625 Inf 1 -2.774 0.1016
## LL Year2022 / NAT Year2023 2.892 0.3181 Inf 1 9.654 <.0001
## LL Year2022 / SIM Year2023 0.691 0.0534 Inf 1 -4.784 <.0001
## NAT Year2022 / SIM Year2022 0.236 0.0280 Inf 1 -12.169 <.0001
## NAT Year2022 / LE Year2023 0.365 0.0461 Inf 1 -7.977 <.0001
## NAT Year2022 / LL Year2023 0.215 0.0256 Inf 1 -12.917 <.0001
## NAT Year2022 / NAT Year2023 0.769 0.1096 Inf 1 -1.842 0.5909
## NAT Year2022 / SIM Year2023 0.184 0.0219 Inf 1 -14.231 <.0001
## SIM Year2022 / LE Year2023 1.545 0.1361 Inf 1 4.935 <.0001
## SIM Year2022 / LL Year2023 0.909 0.0702 Inf 1 -1.234 0.9218
## SIM Year2022 / NAT Year2023 3.259 0.3580 Inf 1 10.753 <.0001
## SIM Year2022 / SIM Year2023 0.779 0.0600 Inf 1 -3.246 0.0258
## LE Year2023 / LL Year2023 0.588 0.0523 Inf 1 -5.971 <.0001
## LE Year2023 / NAT Year2023 2.109 0.2495 Inf 1 6.312 <.0001
## LE Year2023 / SIM Year2023 0.504 0.0447 Inf 1 -7.729 <.0001
## LL Year2023 / NAT Year2023 3.585 0.3958 Inf 1 11.563 <.0001
## LL Year2023 / SIM Year2023 0.857 0.0667 Inf 1 -1.988 0.4901
## NAT Year2023 / SIM Year2023 0.239 0.0264 Inf 1 -12.981 <.0001
##
## P value adjustment: tukey method for comparing a family of 8 estimates
## Tests are performed on the log odds ratio scale
##
## Call:
## glm(formula = cbind(seeded_cover, tot_cover) ~ Treatment + Year,
## family = binomial, data = Total_Cover_Early)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.8127 -1.5620 -0.8626 0.7913 5.1753
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.9526 0.1056 -27.955 < 2e-16 ***
## TreatmentLL -3.3561 0.3411 -9.838 < 2e-16 ***
## TreatmentNAT -2.0984 0.1995 -10.519 < 2e-16 ***
## TreatmentSIM -0.3623 0.1070 -3.387 0.000707 ***
## Year2023 1.0821 0.1106 9.780 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 734.86 on 47 degrees of freedom
## Residual deviance: 247.79 on 43 degrees of freedom
## AIC: 381.05
##
## Number of Fisher Scoring iterations: 6
## Analysis of Deviance Table (Type III tests)
##
## Response: cbind(seeded_cover, tot_cover)
## LR Chisq Df Pr(>Chisq)
## Treatment 373.50 3 < 2.2e-16 ***
## Year 106.69 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## contrast odds.ratio SE df null z.ratio p.value
## LE / LL 28.6759 9.7820 Inf 1 9.838 <.0001
## LE / NAT 8.1535 1.6265 Inf 1 10.519 <.0001
## LE / SIM 1.4366 0.1537 Inf 1 3.387 0.0039
## LL / NAT 0.2843 0.1089 Inf 1 -3.284 0.0056
## LL / SIM 0.0501 0.0172 Inf 1 -8.706 <.0001
## NAT / SIM 0.1762 0.0360 Inf 1 -8.502 <.0001
##
## Results are averaged over the levels of: Year
## P value adjustment: tukey method for comparing a family of 4 estimates
## Tests are performed on the log odds ratio scale
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(seeded_cover, tot_cover) ~ Treatment + Year + (1 | Block)
## Data: Total_Cover_Late
##
## AIC BIC logLik deviance df.resid
## 605.3 616.5 -296.6 593.3 42
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.3147 -1.7714 -0.6057 1.2574 6.5425
##
## Random effects:
## Groups Name Variance Std.Dev.
## Block (Intercept) 0.008129 0.09016
## Number of obs: 48, groups: Block, 6
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.74863 0.09364 -29.352 < 2e-16 ***
## TreatmentLL 1.59358 0.09078 17.554 < 2e-16 ***
## TreatmentNAT 0.16919 0.11161 1.516 0.12953
## TreatmentSIM 1.55200 0.09191 16.885 < 2e-16 ***
## Year2023 0.13908 0.05056 2.751 0.00594 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TrtmLL TrtNAT TrtSIM
## TreatmentLL -0.796
## TreatmntNAT -0.643 0.664
## TreatmntSIM -0.789 0.807 0.656
## Year2023 -0.279 0.018 -0.002 0.024
## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: cbind(seeded_cover, tot_cover)
## Chisq Df Pr(>Chisq)
## (Intercept) 861.5598 1 < 2e-16 ***
## Treatment 566.5093 3 < 2e-16 ***
## Year 7.5683 1 0.00594 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## contrast odds.ratio SE df null z.ratio p.value
## LE / LL 0.203 0.0184 Inf 1 -17.554 <.0001
## LE / NAT 0.844 0.0942 Inf 1 -1.516 0.4278
## LE / SIM 0.212 0.0195 Inf 1 -16.885 <.0001
## LL / NAT 4.155 0.3536 Inf 1 16.738 <.0001
## LL / SIM 1.042 0.0592 Inf 1 0.732 0.8841
## NAT / SIM 0.251 0.0216 Inf 1 -16.031 <.0001
##
## Results are averaged over the levels of: Year
## P value adjustment: tukey method for comparing a family of 4 estimates
## Tests are performed on the log odds ratio scale
PERMANOVA analysis revealed that treatment did have a significant effect on plant community composition. Treatment also explained ~ 20% of the variation in the plant community. Post-hoc analysis shows that LE, NON, and NAT treatments were not significantly different from each other but did differ in plant community composition compared to the SIM and LL seeding treatments. Additionally, the SIM and LL treatments were only marginally different (p = 0.053). Based on the NMDS above, the LL plant community appears as a subset of the SIM plant community. Despite receiving the same pool of seeded species, species added to the plots later in the growing season did not establish as well. Barriers to establishment also appeared higher for spring/summer dispersing species compared to fall dispersing species, particularly when added without priority.
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 1000
##
## adonis2(formula = inner.max ~ Treatment + Year, data = inner_wide, permutations = 1000, method = "bray")
## Df SumOfSqs R2 F Pr(>F)
## Treatment 4 3.1762 0.17687 3.3089 0.000999 ***
## Year 1 1.8224 0.10148 7.5940 0.000999 ***
## Residual 54 12.9588 0.72164
## Total 59 17.9574 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## $parent_call
## [1] "inner.max ~ Treatment + Year , strata = Null , permutations 999"
##
## $LE_vs_LL
## Df SumOfSqs R2 F Pr(>F)
## Treatment 1 0.9839 0.14313 4.2642 0.001 ***
## Year 1 1.0450 0.15202 4.5291 0.001 ***
## Residual 21 4.8456 0.70486
## Total 23 6.8745 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $LE_vs_NAT
## Df SumOfSqs R2 F Pr(>F)
## Treatment 1 0.4410 0.06307 1.6582 0.040 *
## Year 1 0.9656 0.13812 3.6310 0.001 ***
## Residual 21 5.5845 0.79881
## Total 23 6.9911 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $LE_vs_NON
## Df SumOfSqs R2 F Pr(>F)
## Treatment 1 0.6830 0.10287 2.8438 0.001 ***
## Year 1 0.9126 0.13746 3.8000 0.001 ***
## Residual 21 5.0433 0.75967
## Total 23 6.6389 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $LE_vs_SIM
## Df SumOfSqs R2 F Pr(>F)
## Treatment 1 0.7381 0.10692 2.9548 0.001 ***
## Year 1 0.9197 0.13323 3.6819 0.001 ***
## Residual 21 5.2454 0.75986
## Total 23 6.9031 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $LL_vs_NAT
## Df SumOfSqs R2 F Pr(>F)
## Treatment 1 0.7043 0.10632 2.8839 0.001 ***
## Year 1 0.7917 0.11951 3.2417 0.001 ***
## Residual 21 5.1287 0.77418
## Total 23 6.6247 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $LL_vs_NON
## Df SumOfSqs R2 F Pr(>F)
## Treatment 1 1.1298 0.17501 5.2743 0.001 ***
## Year 1 0.8276 0.12820 3.8636 0.001 ***
## Residual 21 4.4985 0.69680
## Total 23 6.4560 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $LL_vs_SIM
## Df SumOfSqs R2 F Pr(>F)
## Treatment 1 0.5752 0.09413 2.6448 0.002 **
## Year 1 0.9681 0.15842 4.4511 0.001 ***
## Residual 21 4.5672 0.74744
## Total 23 6.1105 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $NAT_vs_NON
## Df SumOfSqs R2 F Pr(>F)
## Treatment 1 0.4520 0.07021 1.8176 0.033 *
## Year 1 0.7635 0.11859 3.0701 0.001 ***
## Residual 21 5.2222 0.81119
## Total 23 6.4377 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $NAT_vs_SIM
## Df SumOfSqs R2 F Pr(>F)
## Treatment 1 0.9454 0.13241 3.6145 0.001 ***
## Year 1 0.7019 0.09831 2.6836 0.002 **
## Residual 21 5.4929 0.76928
## Total 23 7.1403 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $NON_vs_SIM
## Df SumOfSqs R2 F Pr(>F)
## Treatment 1 1.2878 0.18695 5.5308 0.001 ***
## Year 1 0.7109 0.10321 3.0533 0.001 ***
## Residual 21 4.8897 0.70984
## Total 23 6.8884 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## attr(,"class")
## [1] "pwadstrata" "list"
# Precipitation
month_tot <- weather %>%
filter(Year != "2021") %>%
group_by(Month, Year) %>%
summarize(tot_month_PRP_mm = sum(Total_PRP_mm))
## `summarise()` has grouped output by 'Month'. You can override using the
## `.groups` argument.
year_tot <- weather %>%
filter(Year != "2021") %>%
group_by(Year) %>%
summarize(tot_PRP_mm = sum(Total_PRP_mm))
library(mosaic)
## Registered S3 method overwritten by 'mosaic':
## method from
## fortify.SpatialPolygonsDataFrame ggplot2
##
## The 'mosaic' package masks several functions from core packages in order to add
## additional features. The original behavior of these functions should not be affected by this.
##
## Attaching package: 'mosaic'
##
## The following objects are masked from 'package:rstatix':
##
## cor_test, prop_test, t_test
##
## The following object is masked from 'package:lmerTest':
##
## rand
##
## The following object is masked from 'package:lme4':
##
## factorize
##
## The following object is masked from 'package:Matrix':
##
## mean
##
## The following object is masked from 'package:cowplot':
##
## theme_map
##
## The following object is masked from 'package:permute':
##
## shuffle
##
## The following objects are masked from 'package:dplyr':
##
## count, do, tally
##
## The following object is masked from 'package:purrr':
##
## cross
##
## The following object is masked from 'package:ggplot2':
##
## stat
##
## The following objects are masked from 'package:stats':
##
## binom.test, cor, cor.test, cov, fivenum, IQR, median, prop.test,
## quantile, sd, t.test, var
##
## The following objects are masked from 'package:base':
##
## max, mean, min, prod, range, sample, sum
favstats(year_tot$tot_PRP_mm)
## min Q1 median Q3 max mean sd n missing
## 665.988 848.995 934.72 1035.304 1148.842 927.735 143.0993 10 0
# Average air temperature
library(mosaic)
favstats(weather$Average_Air_Temp_C)
## min Q1 median Q3 max mean sd n missing
## -21.5 4.277778 13.33333 21.55556 32.66667 12.39584 10.45545 4018 0